#对ADdata采用log2进行标准化，并利用ANOVA检验分别求所有蛋白在3组之间的P值。

{
  #setwd
  rm(list=ls())
  
  setwd("C:\\Users\\37761\\Desktop\\R语言\\第五次实习")
  
  #read data
  data=read.csv("C:\\Users\\37761\\Desktop\\R语言\\k.csv")
  rownames(data)=data[,1]
  data=data[,-1]
  
  #group_info
  group_info = colnames(data)
  group_info <- gsub("[0-9]", "", group_info)
  
  #log2
  data1 = data
  data1=log2(data)
  data1[data1==-Inf]=NA
  data1[101,]=group_info
  rownames(data1)[101]='group'  
}

#analysis
{
  #构建结果矩阵
  rows=rownames(data1)#提取gene name list
  p_ending=data.frame(rows)#根据list组成数据框
  
  #下面是循环体
  for (i in 1:100) 
  {
    #提取列数据
    data2=data1[c(i,101),]
    data2=t(data2)
    data2=as.data.frame(data2)
    colnames(data2)[1]='x'
    
    
    #检验单组内na数据是否符合要求
    c1=data2[data2$group=="ad",]
    c2=data2[data2$group=="ctl",]
    c3=data2[data2$group=="asym",]
    
    
    #判断数值(尽量用变量表示,尽量直观书写，看起来比较清晰但其实没用)
    c1c=((39-sum(is.na(c1$x)))<3)
    c2c=((10-sum(is.na(c2$x)))<3)
    c3c=((15-sum(is.na(c3$x)))<3)
    
    
    #判别获取结果(只要存在一组信息出现一组信息非0数据少于三个的时候的，直接赋值NA)
    if(c1c|c2c|c3c)
    {
      p_ending[i,2]=NA#否则进运算
    }
    if((c1c|c2c|c3c)==F)#这样写直观点
    {
      data2$x=as.numeric(data2$x)
      fit<-oneway.test(x~group_info,data=data2)
      p_ending[i,2]=fit[["p.value"]]
    }
    
  }
  
  #整理结果
  p_ending=p_ending[-1,]
  colnames(p_ending)=c('gene','p')
  
}

#show 
p_ending


######GO富集分析，参考万舸帆

rm(list=ls())
load("C:\\Users\\37761\\Desktop\\R语言\\第四次\\volcano.RData")

table=prostat[prostat$P<0.05,]

#analysis(无特殊要求，固不转换为entrez id处理，选择keytype sumbol进行分析)
ego<- enrichGO(table$ID,OrgDb="org.Hs.eg.db",keyType = "SYMBOL",ont = "all", qvalueCutoff = 0.05,readable = TRUE) 
###设置背景geneset时，结果不理想，可视化失败固不进行额外i展示

#plot(无特殊要求，固进行简单可视化，按照默认设置参数进行绘图)
dotplot(ego,showCategory=10)

png(file="dot_ego.png",width = 1000,height = 1000)
dotplot(ego,showCategory=10)
dev.off()

